AD-A139  552  SOLUTION  TECHNIQUES  IN  FINITE  ELEMENT  ANALVSIS(U) 
CALIFORNIA  UNIV  BERKELEV  DEPT  OF  CIVIL  ENGINEERING 
B  NOUR-ONID  ET  AL.  NAV  83  NCEL-CR-82.  027 
UNCLASSIFIED  N62474-82-C-8277  F/G  20/11 


MICROCOPY  resolution  test  chart 

huiqml  wunru  or  tuxorms 


MIC  FILE  COPY  AD  A1  395  52 


CR  83.027 


NAVAL  CIVIL  ENGINEERING  LABORATORY 
Port  Hueneme,  California 


Sponsored  by 

NAVAL  FACILITIES  ENGINEERING  COMMAND 


SOLUTION  TECHNIQUES  IN  FINITE  ELEMENT  ANALYSIS 


May  1983 


An  Investigation  Conducted  by 
UNIVERSITY  OF  CALIFORNIA 
Department  of  Civil  Engineering 
Berkeley,  California 


N62474-82-C-8277 


I 

I 

1 

f 


Approved  for  public  release-,  distribution  unlimited 


METRIC  CONVERSION  FACTORS 


t  in  -  2.S4  («Jt*ct'v)  For  Olhtfr  «Ksct  conv*r*>ont  »nd  mor*  detailed  tet>t«.  W«  NBS 
Muc  PuW  286,  Unit*  o*  Weights  and  Maesuret,  Price  S2.25.  SO  Cetalofl  No.  Cl  3.10.286. 


Unclassi  fied 


siCu«itv  classification  o 


REPORT  DOCUMENTATION  PAGE 


REPORT  NUMSFR 

CR  83.027 


READ  INSTRUCTIONS 
BEFORE  COMPLETING  FORM 


i  RECIPIENT’S  CATALOG  NUMBER 


4  TiTlC  fmnd  Submit)  S  TYPE  QF,  REPORT  A  PERIOD  COVER 

TECHNIQUES  IN  FINITE  ELEMENT  Sfr  1982  -  Jan  1983 


S  TYPE  OF  REPORT  A  PERIOD  COVERED 


ANALYSIS 


6  PERFORMING  ORG  REPORT  NUMBER 


•  CONTRACT  OR  GRANT  NUMBERrO 


Bahrain  Nour-Omid,  Carlos  Rodrigues,  and  do  r  a  on 

Robert  L.  Taylor  N62474-82-C-8277 


*  PERFORMING  ORGANIZATION  NAME  ANO  AODRESS 

UNIVERSITY  OF  CALIFORNIA 
Department  of  Civil  Engineering 
Berkeley,  Calif. 


It  CONTROLLING  OFFICE  NAME  AND  ADDRESS 

Naval  Civil  Engineering  Laboratory 
Port  Hueneme,  Calif. 


to  PROGRAM  ELEMENT.  PROJECT.  TASK 
AREA  A  WORK  UNIT  NUMBERS 

YR02  3.03.01. 003 


>2  REPORT  DATE 

May  1983 


U  NUMBER  OF  PAOES 


TJ  MONITORING  AGENCY  NAME  A  AOORESSfi/  dilltrant  from  Controlling  Olhct)  IS  SECURITY  CLASS  (ot  t hit  rtport) 

Naval  Facilities  Engineering  Command 

?nn  <:+nUaii  ci.roo*  Unclassi  fied _ 

C\J\J  JlUVa  I  I  Jtrecl  1S«  DECLASSIFICATION  DOWNGRADING 

Alexandria,  Va.  sc"tDuLt 


6  DISTRIBUTION  STATEMENT  folthi*  Rtporr, 


Approved  for  public  release;  distribution  unlimited. 


17  DISTRIBUTION  STATEMENT  (of  tht  tbatrtci  anttrad  in  Block  20.  >1  dilltrant  from  Rtport) 


19  FEY  WOROS  ( Conrinu*  on  rtvtrtt  tidt  if  ntcttttry  and  tdanify  6v  block  nutnbar) 

Structural  mechanics;  finite  elements;  nonlinear  algebraic 
equations;  numerical  solution  methods 


20  ABSTRACT  (Conttnut  on  ra»atma  ••  da  If  nacaaaary  and  tdannfy  by  Alocfc  nittnbar) 


,  A  desirable  advantage  of  iterative  methods  is  that  they 
provide  means  of  controlling  the  accuracy  of  the  solution.  In 
particular,  when  low  levels  of  accuracy  are  required  this  can 
result  in  faster  algorithms  than  the  direct  methods.  The  use  of 
the  conjugate  gradient  algorithm  to  solve  the  linearized  system 
of  equations  is  considered.  A  preconditioning  matrix  based  on  a 


,  1473  EDITION  OF  1  NOV  *5  IS  OBSOLETE 


Unclassi fied _ 

SECURITY  CLASSIFICATION  OF  TmiS  PAGE  D*f*  Enfer*<#) 


Unc I  ass i fied 


SCCwAifv  Cc  ASWiC  At«0*  ' 


T*»‘V  P*r4|  r »»«•«•  D#i •  Pm'0,,u 


^splitting  method  is  constructed.  The  outcome  is  an  algorithm 
wnch  results  in  substantial  reduction  in  storaqe  over  direct 
methods.  The  above  method  is  compared  with  its  rivals  on  several 
quite  different  problems  in  structural  mechanics  and  favorable 
results  were  obtained. 


FORM 

»  JAM*  » 


CO'TiO*  or  *  MOV  ft*  (ft  0»*OLtTt 


_ Unclassi  fied _ 

11  CU*|Tv  Cl»»Si'iCAT,ON  Of  T«|J  P*ol  (»».n  c«r» 


TABLE  OF  CONTENTS 


1.  Introduction  . 

2.  A  Preconditioned  Conjugate  Gradient  Method 

3.  Splitting  Methods  . 

4.  Solution  of  Static  Problems  . 

5.  Solution  of  Dynamic  Problems  . 

6.  Numerical  Examples  . 

6. 1  Static  Examples:  Linear  Elastic  . 

6.2  Static  Examples:  Nonlinear  Elastic  . 

6.3  Dynamic  Examples  . 

7.  Closure  . 

References  . 


WHS 

dtic  t:s 

U?w>naoi;nce4 

Justificatlw; 


Distribution/ 

Availability 

and/or 

Dist  '  Special 

Al  f 


auD 


1.  Introduction 


The  finite  element  method  of  discretization  is  used  to  reduce  many  complex 
continuum  problems  to  discrete  systems.  Although  this  reduction  is  the  most 
important  step  in  the  overall  analysis  of  a  structure,  solving  the  discrete  prob¬ 
lem  is  far  from  trivial.  In  general,  the  reduced  system  is  nonlinear  and  an  itera¬ 
tive  method  must  be  employed  to  arrive  at  the  solution.  Most  solution  methods 
are  based  on  some  form  of  Newton's  method  in  which  the  nonlinear  problem  is 
linearized  using  an  initial  approximation  to  arrive  at  a  linear  set  of  simultaneous 
equations.  The  solution  of  the  set  of  linear  equations  leads  to  a  correction  of  the 
initial  approximation.  When  solving  the  linear  equations,  one  should  not  loose 
sight  of  the  primary  objective:  solving  the  nonlinear  problem. 

Iterative  methods,  such  as  the  conjugate  gradient  or  Lanczos  method,  are 
among  the  many  methods  that  may  be  used  to  solve  systems  of  linear  equations. 
The  advantage  of  these  methods,  when  used  as  the  inner  loop  of  the  Newton 
iteration,  is  twofold. 

(i)  The  linear  equation  may  be  solved  to  any  desired  level  of  accuracy  as 
governed  by  the  Newton  iteration. 

(ii)  A  considerable  reduction  in  storage  can  be  achieved  when  no  triangular  fac¬ 
torization  need  be  performed. 

In  [4]  a  method  was  developed,  based  on  the  preconditioned  Lanczos 
method,  to  realize  some  of  the  advantages  of  iterative  methods.  In  this  previous 
study,  the  triangular  factors  of  the  initial  tangent  matrix  were  used  to  form  a 
preconditioning  matrix  for  the  subsequent  solution  steps.  In  the  present  we 
have  eliminated  factorizations  by  employing  other  preconditioners  and  further, 
have  reduced  the  storage  needs  of  the  method. 
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2.  A  Preconditioned  Conjugate  Gradient  Method 

An  essential  step  in  nonlinear  analysis  of  structures  is  solving  a  linear  sys¬ 
tem  of  algebraic  equations.  The  preconditioned  conjugate  gradient  method 
(hereafter  called  PCG)  is  one  of  the  many  procedures  for  solving 

r  =  b  -  Ax  =  0  (2.1) 

where  A  is  an  nxn  symmetric  positive  definite  matrix  (which  is  sparsely  popu¬ 
lated)  and  b  is  the  right-hand  side  vector.  In  the  case  of  static  analysis,  A  is  the 
current  tangent  matrix  and  in  the  case  of  dynamic  analysis,  A  depends  on  the 
mass,  damping  and  stiffness  matrices,  as  well  as  the  time  increment. 

The  initial  popularity  of  the  conjugate  gradient  method  was  due  to  a 
number  of  factors.  In  exact  arithmetic  the  method  required  a  maximum  of  n 
iterations  to  solve  (2.1)  which  made  the  method  superior  to  other  iterative 
methods,  in  fact  conjugate  gradient  is  in  the  class  of  semi-iterative  methods 
which  also  includes  the  Lanczos  algorithm  [6].  The  disadvantage  of  direct 
methods  is  their  large  storage  demands  for  keeping  the  factors  of  A  The  only 
interface  between  the  conjugate  gradient  method  and  A  is  through  the  product 
Av  for  a  given  vector  ▼.  This  is  an  elegant  way  of  taking  advantage  of  sparsity  of 
A  which  has  the  added  advantage  that  A  need  not  be  known  explicitly  but  only  a 
means  of  computing  the  matrix-vector  product  is  required. 

The  popularity  of  the  conjugate  gradient  method  vanished  once  it  was  found 
that  under  certain  conditions  the  method  required  as  many  as  Sn  or  6n  steps  to 
reduce  the  residual  to  the  desired  level.  This  blemish  is  accounted  for  by  the 
strong  influence  of  round-off  error. 

The  addition  of  preconditioning  eliminated  this  difficulty.  Instead  of  solving 
(2.1)  we  solve 


P'Ax  =  P'b 


(22) 
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for  some  appropriate  choice  of  P.  The  object  then  is  to  choose  P  such  that  the 
coefficient  matrix  of  (2.2)  is  well  conditioned. 

Theoretical  considerations  suggest  that  at  the  end  of  each  iteration  of  CG 

V7  „  i 

the  residual  norm  is  reduced  by  a  factor  ■-  ■  when  solving  (2.1)  where  k  is 

the  condition  number  of  A.  defined  by  k  =  ||A||  {| A- 1 II -  See  [l]  for  more  details. 
Note  that  when  *  =  1.  one  iteration  is  sufficient  to  solve  the  equation.  This  pro¬ 
vides  us  with  a  guideline  for  choosing  P.  For  a  well  chosen  P  only  a  few  iterations 
reduce  the  residual  norm  to  the  desired  level.  Here  we  give  an  outline  for  the 
preconditioned  conjugate  algorithm: 

Given  an  initial  guess  Xo.  a  positive  definite  preconditioning  matrix  P,  the  matrix 
A  and  the  right  hand  side  be 

(1)  Set  pt,  =  r0  =  b  -  Axe 

(2)  Solve  Pdo  =  r0.  for  do 

(3)  for  k  -  0.  1.  2.  •  •  •  until  convergence  do 

(a)  a*  =  (r*<d*)/(lfe.Apfc) 

(b)  it  *  i  =  *  +  a*pk 

(c)  r*tl  =  rk-a*Ap* 

(d)  Solve  Fdn-!  =  r*»i 

(e)  flu  *  (r*  +  |. d*+i)/  (r*,dk ) 

(0  =  4s>i+^*Pfc 

The  operation  (e.u)  denotes  the  inner  product  Tru.  The  algorithm  gen¬ 
erates  a  sequence  of  approximations  to  the  solution  x  with  a  corresponding  resi¬ 
dual  vector  r*.  The  termination  criterion  can  be  chosen  based  on  these  quanti¬ 
ties.  In  addition  to  storage  demands  for  A  and  P  the  algorithm  requires  storage 
for  4  vectors.  The  total  number  of  operation  per  iteration  is  NZA  +  NZP  +  bN . 
where  NZA  and  NZM  are  the  number  of  operations  for  forming  Au  and  P~'v  for  a 
given  n  end  ▼. 


3.  Splitting  Methods 


Here  we  turn  to  a  topic  which  at  first  sight  may  seem  unrelated  to  the  solu¬ 
tion  of  nonlinear  algebraic  equations.  Consider  the  system  of  first  order 
differential  equations 

i  =  f(xf )  (3.1) 

where  x  is  an  n  -dimensional  vector,  the  superposed  dot,  (  ),  denotes 

differentiation  with  respect  to  time  and  f  is  a  function  of  the  unknown  vector  x 

and  t. 

We  consider  a  special  form  of  f  which  can  be  written  as  a  sum  of  its  subcom¬ 
ponents  f*. 

*=  tu  (3.2) 

t  -  l 

Under  these  conditions  the  original  problem  can  be  thought  &s  a  sum  of  s  sub- 
problems 

X  =  (<(&*)  i  =  l . s  (3.3) 

In  the  case  of  finite  element  discretization  of  the  spatial  domain  the  sum  in  (3.2) 

ranges  over  the  elements.  In  other  cases  the  splitting  may  be  formed  by  other 

means,  one  of  which  is  demonstrated  in  the  following  section. 

A  consistent  algorithm  for  the  solution  of  (3.1).  based  on  the  notion  of  a 
splitting  technique  [3],  can  now  be  constructed  as  a  product  of  algorithms  for 
the  sub-problems.  In  other  words,  write  the  algorithm  for  (3.3)  as 

=  S^ix,*]  (3.4) 

where  S/**  is  an  operator  denoting  the  algorithm  and  the  index  m  ranges  over 

the  increment  in  time.  h.  Then  the  algorithm  for  (3.1)  can  be  written  as 


(3.5) 

=  f[ 

(36) 

i  ■  1 
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One  of  the  disadvantages  of  the  splitting  method  is  its  low  accuracy.  The 
best  that  these  methods  can  achieve  is  second  order  accuracy.  That  is  the  trun¬ 
cation  error  is  of  the  order  of  h*  at  best.  In  the  sequel  we  will  use  the  above  pro¬ 
cedure  to  construct  a  preconditioning  matrix  for  the  conjugate  gradient  algo¬ 
rithm  described  in  section  2.  The  inherent  inaccuracy  of  the  splitting  method 


poses  no  problem  since  the  algorithm  is  used  only  as  a  preconditioner  and 
therefore  one  can  obtain  very  high  accuracies  through  the  conjugate  gradient 
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4.  Solution  of  Static  Problems 

Consider  the  system  of  linear  first  order  differential  equations 

tx  +  Ax  —  b  (4.1) 

where  r  is  a  given  parameter.  Formally  the  solution  to  equation  (4.1)  is 

tit)  =  e-*/T(xo  -  A-1b)  +  A-,b  (4.2) 

where  Xo  =  x(0).  is  an  initial  condition.  We  observe  from  (4.2)  thet  as  t  ends  to 

infinity  x(f )  converges  to  the  solution  of  (2.1)  for  r  >  0.  Consiquently  (4.1)  may 
be  utilised  to  solve  the  linear  equations  (2.1).  Indeed  this  approach  has  been 
suggested  previously  (e.g.,  see  [5]).  In  general  the  exponential  of  a  large  matrix 
cannot  be  easily  computed  and  a  numerical  solution  of  (4.1)  must  be  used.  In 
order  to  achieve  a  soluion  of  (2.1)  a  numerical  solution  to  (4.1)  must  be  assymp- 
totically  correct  for  infinite  Af,  or  a  very  large  number  of  time  steps  must  be 
used  to  compute  the  solution  at  infinite  time.  Here  we  are  not  concerned  with 
constructing  accurate  solution  to  (4.1)  rather  we  consider  the  method  as  a 
means  of  constructing  a  suitable  preconditioning  matrix  for  the  conjugate  gra¬ 
dient  algorithm  described  above. 

Splitting  methods  may  be  applied  to  any  problem  of  the  form 

i=Bx  (4.3) 

where  B  is  an  additive  operator  defined  by 

B  =  £  B<  (4.4) 

t  *  > 

such  that  the  equations 

x  =  BtX  i  =  1 . s  (4.5) 

are  significantly  easier  to  solve  than  the  original  equations.  The  time  stepping 

algorithm  for  the  global  problem  is  then  the  product  of  ail  the  time  stepping 

algorithms  for  the  subproblems  with  a  fractional  time  step  h/ s  [3]. 

The  coefficient  matrix  A  in  (2.1)  may  be  written  as  the  sum  of  its  diagonal 
matrix,  D,  a  strictly  lower  triangular  matrix,  L  and  a  srictly  upper  triangular 
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matrix  such  that 

A  =  0$D  +  L)  +  (&D  +  L)r  (4.6) 

The  associated  subproblems,  x  =  -0$D  +  L)x  and  x  =  -0$D  +  Lr)x  can  be  solved 

easily.  Applying  a  backward  difference  method  with  a  time  step  h/Z  to  each  of 

the  subproblems,  we  arrive  at 

|-i|  j-ir 

**♦1  =  l*  2l  I+£{D+2Lr)  £b+x«  (4.7) 

where  is  an  approximation  to  x(mA).  For  initial  condition  Xo  =  0  we  get  an 
approximation  to  x(A) 

*i  =  {^rjl  +  £<D  +  2L)  I  +  £(D  +  21/)  b  (4.8) 

which  is  compared  to  the  exact  solution 

x(M  =  [i  -  e"i‘/T)A-1  b  (4.9) 

Comparing  equations  (4.8)  and  (4.9)  suggests  that  the  coefficient  matrix  in  (4.8) 

may  be  a  good  approximation  to  A-1  for  large  h.  and  may  therefore  be  an 
effective  preconditioning  matrix.  The  scalar  factor  may  be  ignored  for 
preconditioning  purposes. 

When  using  this  in  conjunction  with  conjugate  gradient  algorithm  of  section 
2  the  preconditioning  matrix  becomes 

P  =  (1  +  «/  2D  +  oL)(I  +  u/  2D  +  t>Lr)  (4. 10) 

where  u  =  h/ 2t  is  now  a  free  parameter. 

To  simplify  the  choice  of  u  we  scale  the  stiffness  matrix  A  such  that  diagonal 
of  A  is  unity.  The  resulting  matrix  is  A  =  D^AD"#.  The  system  of  equations  (2.1) 
now  becomes 

At  =  b  (4.11) 

where  t  =  D*x  and  b  =  D~*b. 

The  preconditioned  matrix  must  now  be  applied  to  (4.11)  resulting  in 


P  =  (I  ♦  6>L)(I  +  wtT)  (4.12) 

where  A  =  L  +  Lr.  It  is  easy  to  show  that  preconditioning  (4.11)  using  P  is 


equivalent  to  preconditioning  (2.1)  with 

P  =  (D  +  oL)D-‘(D+  «Lr)  (4.13) 

This  can  be  identified  as  a  member  of  the  class  of  incomplete  Choleski  precondi¬ 
tioners  [2).  Note  that  when  u  =  0,  P  becomes  the  diagonal  matrix  D,  resulting  in 
the  simplest  form  of  preconditioning;  diagonal  scaling.  When  u  =  1  then 
P  =  A  +  LD-1Lr.  The  error  matrix  LD_,Lr  is  rank  deficient  since  L  has  zero  diag¬ 
onals.  If  the  norm  of  D  is  larger  then  the  norm  of  L  then  the  norm  of  the  error 
m-lrix  will  be  small  compared  to  the  norm  of  A.  consequently,  for  most  prob¬ 
lems  it  is  expected  that  the  optimum  u  will  be  close  to  unity. 
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5.  Solution  of  Dynamic  Problems 

We  construct  a  preconditioning  matrix  for  the  linear  system  of  equations 
arising  in  a  step-by-step  algorithm  for  dynamic  analysis  of  linear  and  nonlinear 
structures.  In  particular,  we  consider  the  Newmark  algorithm  and  the  precondi¬ 
tioning  matrix  follows  from  the  splitting  method  of  section  3,  in  much  the  same 
way  as  for  the  static  problem. 

Consider  the  linear  equations  of  motion 

Mii  +  Ku=r  (5.1) 

where  M  is  the  diagonal  mass  matrix,  K  is  the  stiffness  matrix,  f  is  the  external 

load  vector  and  u  is  the  response  of  the  structure.  For  simplicity,  we  ignore  the 

damping  effects,  but  all  the  following  results  may  be  extended  easily  to  the 

damped  case.  The  linear  system  of  equations  arising  at  every  time  step  of  the 

Newmark  method  is 


Ax  =  b  (5.2) 

where 


and 


A  =  K  + 


(5.3) 


b  =  ft  ♦  ♦  ~^-M[u,  +  ***  +  OS  -  PW 1  (5.4) 

Here  ▼  and  a  are  velocity  and  acceleration  vectors,  respectively,  Af  is  the 
specified  time  increment,  t  is  the  time  and  x  is  the  increment  of  displacement 
response.  The  Newmark  parameters  are  chosen  such  that  0  at  ()$  +  y)8/  4  with 
7  as  )$.  The  discretization  in  time  are 

u<  =  «i  ♦  A/W|+  ♦HAf8[(l-2/J)a<  ♦  20a,  )]  (5.5) 

▼,  ♦  At  =  *t  ♦  Af  (1  -  7)a,  ♦  jLt  a,  ♦  tt  (5  8) 

The  object  is  to  solve  (5.2)  without  forming  the  factors  of  A. 

A  splitting  method  similar  to  the  one  used  for  equation  (4.1)  can  now  be 
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applied  to  equation  (5.1).  The  matrix  resulting  from  the  splitting  algorithm  can 
then  be  used  as  a  preconditioner  for  (5.2).  Consider 


p  =  <L  +  ^ 

where  K  =  L  +  Lr.  Multiplying  out  the  terms  in  (5.7).  we  obtain 


=  +  A] 

=  -t-mu*)  +  A]  (5.8) 

where  E(Ata)  =  0Atzm_,Lr. 


The  preconditioned  conjugate  gradient  algorithm  of  section  2  is  invariant 
under  the  scaling  of  the  preconditioning  matrix.  Therefore,  (5.B)  shows  that  P 
will  tend  quadratically  to  the  dynamic  stiffness  matrix  A  as  the  time  step  dimin¬ 
ishes.  In  other  words.  E  tends  to  the  zero  matrix  quadratically  in  At.  We  see 
later  that  this  characteristic  results  in  an  effective  preconditioning  and  the  solu¬ 
tion  of  equation  (5.2)  is  obtained  in  as  few  as  2  or  3  iterations  of  the  precondi¬ 
tioned  conjugate  gradient  algorithm  with  moderately  small  At. 


».v  rr»r  n  ■  nj, ■■  m  ».ii  ■i,i  »y  >i  »;i»»  ■  ■  yi 


.V.  TT-  r.'F. 
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6.  Numerical  Examples 

In  the  following,  we  present  a  few  numerical  examples  to  illustrate  some  of 
the  characteristics  of  the  proposed  preconditioning  matrices  for  the  PCG  algo¬ 
rithm.  This  algorithm  is  implemented  in  FEAP,  a  Finite  Element  Analysis  Pro¬ 
gram  (see  [8],  chapter  24  for  more  details).  All  the  numerical  tests  were  carried 
out  on  a  VAX  11/780  at  the  University  of  California,  Berkeley,  using  double  preci¬ 
sion  computation. 


We  first  present  the  results  to  some  static  analyses,  both  linear  and  non¬ 
linear.  In  these  examples  we  choose  a  stopping  criterion  based  on  the  residual 
vector  and  the  algorithm  is  terminated  as  soon  as  the  norm  of  this  vector  is 
reduced  by  a  factor  smaller  than  a  specified  tolerance.  In  our  calculations  we 
set  the  tolerance  to  10-8.  Next  we  demonstrate  our  algorithm  on  a  few  dynamic 
problems.  The  termination  criterion  is  similar  to  the  static  case  with  a  range  of 
different  tolerances  to  demonstrate  the  effectiveness  of  the  algorithm. 

6.1  Static  Examples:  Linear  elastic 
a)  132  degree-o /-freedom  building 

The  object  of  the  first  problem  is  to  determine  the  influence  of  the  precon¬ 
ditioning  parameter,  <y,  in  eq.  (4.13).  The  total  number  of  PCG  iterations 
required  to  achieve  convergence,  varies  considerably  with  u.  To  illustrate  this 
dependence,  we  chose  the  example  model  shown  in  figure  1  which  is  a  132 
degrees  of  freedom,  multistory  building,  discretized  by  176  2-node  truss  ele¬ 
ments  each  with  the  same  Young's  modulus  (30xl0a).  The  cross-sectional  area 
of  the  girders,  columns  and  diagonals  are  20,  40  and  1  respectively.  A  single  load 
at  the  top  is  applied,  as  shown  in  figure  7. 

In  figure  2,  we  indicate  the  number  of  PCG  iterations  needed  to  converge  as 
a  function  of  the  preconditioning  parameter.  The  shape  of  this  curve  is  charac¬ 
teristic  of  the  proposed  PCG  algorithm  and  consists  of  three  zones: 
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1)  Small  u:  The  preconditioning  matrix  approaches  the  diagonal  matrix  D-  In 
this  case,  the  total  number  of  iterations  is  less  than  that  for  diagonal 
preconditioning. 

2)  Optimum  u:  With  this  value  the  algorithm  takes  the  least  number  of  itera¬ 
tions  to  obtain  the  solution.  Note  that  the  curve  is  quite  flat  around  a>op< 
and  therefore  the  total  number  of  PCG  steps  is  insensitive  to  small  changes 
in  the  value  of  cj.  Further,  as  predicted  before,  the  optimum  u  is  close  to 
unity. 

3)  large  u:  In  this  range  the  preconditioning  matrix  approaches  UI_,Lr  which 
is  a  singular  matrix  (diagonals  of  L  are  zero).  In  this  example  with  u  >  3.0 
the  solution  may  loose  accuracy  in  all  significant  digits  and  eventually 
floating-point  overflow  occur. 

Figure  3  shows  the  evolution  of  the  residual  norm,  ||  rt  || .  normalized  versus 
||r0||.  at  the  t-th  iteration  of  the  PCG  algorithm.  The  residual  at  each  iterate 
exhibits  characteristics  typical  of  conjugate  gradient  method.  Namely,  residual 
norm  remains  large  for  a  relatively  large  number  of  steps  before  convergence 
occurs  to  the  specified  tolerance.  Part  of  this  behavior  is  due  to  the  loss  of 
orthogonality  among  the  conjugate  vectors.  Poor  preconditioning  can  also  con¬ 
tribute  to  slow  convergence. 

b)  "Cantilever  beam "  type  structures 

From  the  insight  we  have  gained  with  the  preceding  example,  we  now 
proceed  to  answer  the  following  question:  How  to  select  the 

No  easy  analytical  solution  can  be  obtained  to  this  question;  depends 
on  the  spectrum  of  A  which  is  not  known  apriori.  However  an  initial  estimate  of 
unity  as  indicated  in  section  4  is  not  an  unreasonable  choice  for  u.  The  numeri¬ 
cal  test  here  is  to  investigate  the  dependence  of  the  number  of  iterations  on  u. 
An  accurate  upper  bound  to  the  total  number  of  PCG  steps  can  be  obtained  if 


the  condition  number  of  the  preconditioned  matrix,  P  lA.  is  known.  However  the 
condition  number  of  P"‘A  depends  on  u. 

The  examples  we  have  chosen  are  summarized  in  figure  4.  Each  problem  is 
computed  with  a  range  of  c/s  to  obtain  cjBp , .  Figure  5  shows  the  number  of  itera¬ 
tions  as  a  function  of  u,  for  these  examples.  Notice  that  all  the  curves  are 
rather  flat  when  close  to  uopt .  moreover,  that  is  close  to  1.0. 

The  following  table  (1)  gives  the  number  of  iterations  for  both  Uopt  and 

a  -  1.0. 


Table  1.  Comparison  of  the  Number  of  PCG  Iterations  for  Various  c/s. 

The  last  column  of  Table  1  shows  the  ratio  of  total  number  of  iterations  over 
the  number  of  degree-of-freedom.  As  expected  this  ratio  remains  below  unity. 
The  next  to  last  column  shows  the  loss  in  optimality  when  using  u  equal  to  unity. 
Except  for  the  132  degree-of-freedom  system  little  loss  in  computational  effort 
results  from  using  u  equal  to  one. 


6.2  Static  Examples:  Nonlinear  Elastic  Problem 
a)  nonlinear  material  problem 

In  much  the  same  way  as  the  Newton-Lanczos  method  [4],  the  PCG  algo¬ 
rithm  was  implemented  within  a  Newton  loop.  The  resulting  algorithm  possessed 
all  the  properties  of  the  Newton-Lanczos  algorithm  with  the  exception  that  it  is 
restricted  to  positive  definite  matrices.  Simplicity  of  programing  various  altera¬ 
tions  was  the  motivating  factor  in  restricting  attention  to  the  PCG  algorithm. 
Our  primary  objective  is  to  compare  the  PCG  algorithm  with  Newton  and 
modified  Newton  strategies.  For  this  comparison  we  use  the  132  degree-of- 
freedom  truss  building  described  above,  but  modified  to  have  the  same  cross- 
section  for  all  the  members  (A  -  20).  Nonlinearity  is  introduced  by  a  simple 
yield  model  in  the  constitutive  equation. 

A  single  load  is  applied  at  the  top  with  sufficient  magnitude  to  produce  a 
nonlinear  maximum  displacement  of  about  twice  the  maximum  elastic  one.  Fig¬ 
ure  6  shows  the  mesh,  the  deformed  structure  and  the  constitutive  equation 
adopted.  In  Table  II,  we  indicate  the  relative  computational  cost  comparisons 
for  different  methods.  We  modified  all  the  algorithms  mentioned  above  to 
include  a  line  search.  This  was  initially  expected  to  reduce  the  final  cost  of  the 
algorithms;  in  fact,  for  this  problem  the  three  methods  were  more  expensive 
when  a  line  search  was  included. 

Looking  at  the  results  in  Table  II,  it  is  interesting  to  note  that  the  PCG  algo¬ 
rithm  required  only  one  more  nonlinear  step  than  the  Newton  method.  Also,  due 
to  the  fact  that  only  the  nonzero  terms  of  the  stiffness  matrix  are  stored,  the 
coat  for  one  matrix-vector  operation  in  the  PCG  algorithm  is  smaller  than  for 
the  other  methods.  For  this  example  the  number  of  terms  in  the  matrix  stored 
in  profile  form  is  1654,  however  the  PCG  algorithm  requires  only  512  nonzero 
terms.  Therefore  the  cost  of  one  matrix-vector  operation  is  about  a  third  of  a 
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profile  multiply.  More  important  is  the  reduction  in  the  over  all  storage 
demand.  For  this  example  the  storage  is  reduced  to  31%  of  the  amount  required 
for  a  profile  stored  solution. 


Method 

No.  of 
iter. 

No.  of 
LU 

factor. 

No.  of 
function 
evaluations 

— 

No.  of 

matrix-vec.  ! 
operations  { 

Newton 

6 

6 

7 

6 

Mod.  Newton 

215 

1 

216 

215 

PCG 

7 

0 

B 

1B6 

Newton  +  LS 

6 

6 

7 

6 

Mod.  Newton  +  LS 

111 

1 

112 

111 

PCG  +  LS 

7 

0 

B 

176 

Table  11.  Cost  Comparisons  for  different  Nonlinear  Methods  (Truss  Example). 

The  average  number  of  PCG  iterations  was  27.  for  a  preconditioning  param¬ 
eter  u  s  1.5.  This  compares  with  37  PCG  iterations  for  a  linear  problem  with 
much  the  same  structure  (see  section  B.la).  A  lower  tolerance  for  PCG  algo¬ 
rithm  is  used  in  the  earlier  stages  of  the  Newton  loop  which  accounts  for  the 
lower  average  number  of  iterations  (see  [4]  for  more  details). 

In  this  test  the  total  cost  for  the  PCG  algorithm  was  twice  the  Newton  cost. 
However,  this  ratio  is  expected  to  drop  well  below  1.0  for  three-dimensional 
structures  where  the  cost  of  a  factorization  is  large  compared  to  the  matrix- 
vector  operation.  Moreover,  as  noted  previosly  we  require  substantially  less 
storage  space. 

b)  finite  deformation  problem 

In  figure  7,  we  show  a  plane  strain  rubber  block  subjected  to  large  deforma¬ 
tion.  We  employ  a  4-node  element  and  a  Mooney-Rivlin  material  as  described  in 


[7].  The  rubber  block  is  discretized  by  144  elements  (12  x  12  mesh)  with  286 
degrees  of  freedom.  The  stiffness  matrix  stored  in  profile  form  requires  7618 
storage  spaces  of  which  only  1799  are  nonzero;  corresponding  to  a  76%  saving  in 
storage  when  using  PCG. 

The  rubber  block  is  stretched  to  50%  of  its  original  length  in  load  5  steps. 
The  cost  comparison  of  both  the  PCG  and  Newton  algorithms  is  summarized  in 
the  following  table. 


displ. 

u 

Nonlinear  PCG 

Newton 

No.  of 
Iterat. 

Averg.  No. 
Matrix  op. 

No.  of 
Iter. 

0.1 

6 

43 

7 

0.2 

6 

45 

6 

0.3 

7 

47 

6 

0.4 

6 

47 

5 

0.5 

6 

44 

5 

Table  III.  Comparision  of  PCG  Algorithm  and  Newton  Method  (Rubber  Block). 

Again,  both  PCG  and  Newton  require  almost  the  same  number  of  nonlinear 
steps  to  converge.  What  is  more  interesting  is  that  the  number  of  PCG  iterations 
is  quite  constant,  even  for  the  highly  nonlinear  range.  When  comparing  this  test 
with  the  previous  132  degree-of-freedom  building  example,  we  notice  that  the 
number  of  iterations  in  the  PCG  algorithm,  as  expected,  does  not  increase  as 
fast  as  the  number  of  degree-of-freedom. 

6.3  Dynamic  example 


In  this  example,  we  wish  to  indicate  the  effectiveness  of  solving  approxi¬ 
mately,  a  linear  elastic  dynamic  problem  using  the  PCG  algorithm.  This  is  done 
for  a  series  of  time  steps  and  tolerances.  Since  the  PCG  algorithm  involves  no 


factorization  steps  it  can  solve  nonlinear  dynamic  problems  with  the  same 
amount  of  effort  as  for  the  linear  case.  However,  in  this  study  we  select  a  linear 
problem.  To  limit  the  computer  costs,  we  selected  a  structure  having  20  4-node 
plane  stress  elements  defined  in  section  6.1(b).  The  dynamic  problem  consists  of 
releasing  the  structure  from  an  initially  deformed  configuration  and  letting  it 
vibrate  freely.  The  mesh,  material  properties  and  initial  state  are  given  in  figure 
8.  The  time  steps  chosen  are  Af  *  0.5.  0.2,  0.1,  0.05  and  0.025  seconds,  which 
correspond  to  1/  120  <  Af  /  T  <  1/  6.  where  T-  3.0  sec.  is  the  fundamental 
period  of  vibration  of  the  structure.  For  comparison,  the  critical  time  step  for 
an  explicit  analysis  would  be  &t„.  =  0.01  sec.,  for  a  bulk  wave  velocity  of  911/s. 
In  order  to  see  the  effect  of  solving  approximately  the  set  of  equations,  we  use 
the  three  following  tolerances:  tol  =  10-4.  10"a  and  10-1.  In  figure  9.  we  plot  the 
y -displacement  of  node  1  for  the  time  step  A t  =0.1  for  the  three  tolerances;  in 
figure  9,  we  show  the  corresponding  relative  error  e,  i.e. 
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where  6s  is  the  displacement  obtained  using  the  Newmark  method,  6pct  is  the 
corresponding  results  obtained  using  PCG  and  60  =  0.172  is  the  initial  applied 
displacement.  The  results  clearly  show  that  the  tolerance  10_l  is  too  large  and 
leads  to  inaccurate  results.  The  error  using  a  tolerance  of  10~B  is  about  1%,  while 
there  is  no  visible  error  for  fol  =  10-4  (less  than  0.01  percent).  When  we  reduce 
the  time  step  to  At  =  0.05  (half  the  preceding),  the  results  improve  substan¬ 
tially:  while  we  see  no  difference  between  tolerances  10-4  and  10'8,  there  is  only 
1  %  error  when  using  10'1  (fig.  10).  For  smaller  time  steps,  no  differences  are 
seen  in  the  first  five  digits. 


Finally,  figure  11  shows  the  average  number  of  iterations  as  a  function  of 
F/At.  The  reduction  of  the  number  of  iterations  as  At  tends  to  zero  is  quite 
interesting:  for  a  tolerance  of  10~2,  this  number  drops  from  30  to  4  iterations 


This  reduction  is  totally  due  to  the  con- 
o  the  dynamic  stiffness  matrix,  A.  Such 
applications,  e.g.  impact  problems. 
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7.  Closure 

In  this  report  we  have  described  our  initial  efforts  to  construct  a  solution 
method  for  the  algebraic  equations  arising  from  finite  element  solution  of  linear 
and  non-linear  problems.  Both  static  and  dynamic  problems  are  considered. 
For  nonlinear  problems,  Newton's  method  is  used  to  generate  a  sequence  of 
linear  problems.  A  preconditioned  conjugate  gradient  method  is  used  to  solve 
the  linear  set  of  equations.  A  method  for  constructing  an  effective  precondition¬ 
ing  matrix  in  terms  of  an  additive  decomposition  of  the  coefficient  matrix  is 
introduced  separately  for  the  static  and  dynamic  cases.  Several  example  prob¬ 
lems  are  solved  demonstrating  the  features  of  the  proposed  method. 

In  order  to  further  evaluate  the  method  additional  work  is  required.  In  par¬ 
ticular  we  recommend  that  the  conjugate  gradient  part  of  the  algorithm  be 
replaced  by  the  Lanczos  method  as  described  in  [4,6],  This  will  permit  con¬ 
sideration  of  indefinite  problems,  such  as  those  resulting  from  use  of  Lagrange 
multiplier  methods  (e.g., contact  problems,  etc.).  In  addition  it  is  essential  to 
test  the  method  on  larger  problems,  preferably  some  three-dimensional  prob¬ 
lems  where  sparsely  populated  coefficient  matrices  with  rather  large  mean 
column  heights  occur.  Further  analyses  for  significant  non-linear  problems 
should  also  be  performed.  Finally,  some  efforts  to  adaptively  compute  an 
optimal  value  for  u  should  be  explored. 


-  20  - 


REFERENCES 

[1]  Concus,  P..  G.  H.  Golub  and  D.  P.  O’Leary,  "A  Generalized  Conjugate  Gradient 
Method  for  the  Numerical  Solution  of  Elliptic  Partial  Differential  Equations.” 
Sparse  Matrix  Computations,  ed.  by  J.  R.  Bunch  and  D.  J.  Rose,  Academic 
Press,  New  York,  1976. 

[2]  Kershaw,  D.  S.,  "The  Incomplete  Cholesky-Conjugate  Gradient  Method  for 
Solution  of  System  of  Linear  Equations,"  Journal  of  Computational  Physics, 
vol.  28,  pp.  43*65,  1978. 

[3]  Gourlay,  A.  R..  "Splitting  Methods  for  Time  Dependent  Partial  Differential 
Equations,"  The  State  of  the  Art  in  Numerical  Analysis,  ed.  by  D.  Jacobs, 
Academic  Press,  New  York,  1977. 

[4]  Nour-Omid,  B.,  B.  N.  Parlett  and  R.  L.  Taylor,  "A  Newton-  Lanczos  Method  for 
Solution  of  Nonlinear  Finite  Element  Equations,  Computers  and  Structures, 
vol.  16,  No.  1-4,  pp.  241-252,  1962. 

[5]  Park,  K.  C.  and  J.  M.  Housner,  "Semi-Implicit  Transient  Analysis  Procedures 
for  Structural  Dynamic  Analysis,"  International  Journal  for  Numerical 
Methods  in  Engineering ,  vol.  18,  pp.  609-622,  1982. 

[6]  Parlett,  B.  N..  "A  New  Look  at  the  Lanczos  Algorithm  for  Solving  Symmetric 
Systems  of  Linear  Equations.  Linear  Algebraic  Applications,  vol.  29,  pp. 
323-346,  1980. 

[7]  Simo,  J.  C.  and  R.  L.  Taylor,  "Penalty  Formulations  for  Rubber-like  Elasti¬ 
city.  Report  No.  UCB/SESM  81/03,  Department  of  Civil  Engineering, 
University  of  California,  Berkeley,  November  1981. 

[8]  Zienkiewicz,  0.  C.,  77ie  Finite  Element  Method,  Third  Edition  .  Mc-Graw  Hill, 
Inc..  London.  1977. 


MODULE 


10 


E  =30.  106 
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Figure  1.  132  Degrees  of  Freedom  Truss  Building. 


Figure  6.  132  Degrees  of  Freedom  Nonlinear  Truss  Building 
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Figure  11.  Improvement  of  Preconditioning  Matrix  with  Reduction  in  Time  Step. 


